%Genetic algorithm
function output = genalg(input,cost_function)
    bounds = input.bounds;
    n_var = size(bounds,1);                 %number of variables
    n_chr = input.n_chr;                    %number of chromosomes
    r_sel = input.r_sel;                    %selection rate
    n_sur = round(n_chr*r_sel);             %number of survived chromosomes
    r_mut = input.r_mut;                    %mutation rate
    n_mut = round(n_chr*r_mut);             %number of mutations
    r_int = input.r_int;                    %introduction rate
    n_int = round(n_chr*r_int);             %number of introductions
    n_gen = input.n_gen;                    %number of generations
    eps = input.eps;
    
    %First generation
    population = zeros(n_var,n_chr); %preallocation
    for i1 = 1:n_chr
        population(:,i1) = unifrnd(bounds(:,1),bounds(:,2));
    end
    
    cost = zeros(1,n_chr);
    cost_record = zeros(1,n_gen);
    cost(1) =  eps*1.1;
    i1 = 1;
    while i1 <= (n_gen - 1) && cost(1) > eps
        %Cost function calculation
        for i2 = 1:n_chr
            cost(i2) = cost_function(population(:,i2));
        end
        [cost,index] = sort(cost,2);
        cost_record(i1) = cost(1);

        population = population(:,index); %sorting chromosomes

        %Normalizing population values to [0,1] interval
        bounds_low = repmat(bounds(:,1),1,n_chr);
        bounds_high = repmat(bounds(:,2),1,n_chr);
        population = (population - bounds_low)./(bounds_high - bounds_low);

        %Mating
        for i2 = (n_chr - n_sur):n_chr
            r_1 = randi(n_sur); %index of 1st parent
            r_2 = randi(n_sur); %index of 2nd parent
            r_3 = rand; %crossover point
            parent_1 = population(:,r_1);
            parent_2 = population(:,r_2);
            population(:,i2) = r_3*parent_1 + (1 - r_3)*parent_2;
        end

        %Mutating
        for i2 = 1:n_mut
            r_1 = randi(n_chr); %index of mutating chromosome
            r_2 = rand; %crossover point
            mutant = population(:,r_1);
            population(:,r_1) = r_2*mutant + (1 - r_2)*mutant;
        end

        %New chromosomes introduction
        for i2 = 1:n_int
            population(:,i2) = unifrnd(bounds(:,1),bounds(:,2));
        end

        %Denormalizing population values
        bounds_low = repmat(bounds(:,1),1,n_chr);
        bounds_high = repmat(bounds(:,2),1,n_chr);
        population = population.*(bounds_high - bounds_low) + bounds_low;
        
        i1 = i1 + 1;
    end
    
    %Final cost function calculation
    for i2 = 1:n_chr
        cost(i2) = cost_function(population(:,i2));
    end
    [cost,index] = sort(cost,2);
    cost_record(i1) = cost(1);

    population = population(:,index); %sorting chromosomes
    
    output.best = population(:,1);
    output.cost = cost(1);
    output.record = cost_record;
end